Contracting cubeΒΆ
In this demo we simulate a unit cube that is fixed at \(x = 0\) and free at \(x = 1\). We use a transversally isotropic material with fiber oriented in the \(x\)-direction.
[1]:
import dolfin
try:
from dolfin_adjoint import (
Constant,
DirichletBC,
Expression,
UnitCubeMesh,
interpolate,
Mesh
)
except ImportError:
from dolfin import (
Constant,
DirichletBC,
interpolate,
Expression,
UnitCubeMesh,
Mesh
)
[2]:
import pulse
pulse.iterate.logger.setLevel(10)
[3]:
# Create mesh
N = 6
mesh = UnitCubeMesh(N, N, N)
[4]:
# Create subdomains
class Free(dolfin.SubDomain):
def inside(self, x, on_boundary):
return x[0] > (1.0 - dolfin.DOLFIN_EPS) and on_boundary
class Fixed(dolfin.SubDomain):
def inside(self, x, on_boundary):
return x[0] < dolfin.DOLFIN_EPS and on_boundary
[5]:
# Create a facet fuction in order to mark the subdomains
ffun = dolfin.MeshFunction("size_t", mesh, 2)
ffun.set_all(0)
[6]:
# Mark the first subdomain with value 1
fixed = Fixed()
fixed_marker = 1
fixed.mark(ffun, fixed_marker)
[7]:
# Mark the second subdomain with value 2
free = Free()
free_marker = 2
free.mark(ffun, free_marker)
[8]:
# Create a cell function (but we are not using it)
cfun = dolfin.MeshFunction("size_t", mesh, 3)
cfun.set_all(0)
[9]:
# Collect the functions containing the markers
marker_functions = pulse.MarkerFunctions(ffun=ffun, cfun=cfun)
[10]:
# Create mictrotructure
V_f = pulse.QuadratureSpace(mesh, 4)
[11]:
# Fibers
f0 = interpolate(Expression(("1.0", "0.0", "0.0"), degree=1), V_f)
# Sheets
s0 = interpolate(Expression(("0.0", "1.0", "0.0"), degree=1), V_f)
# Fiber-sheet normal
n0 = interpolate(Expression(("0.0", "0.0", "1.0"), degree=1), V_f)
[12]:
# Collect the mictrotructure
microstructure = pulse.Microstructure(f0=f0, s0=s0, n0=n0)
[13]:
# Create the geometry
geometry = pulse.Geometry(
mesh=mesh,
marker_functions=marker_functions,
microstructure=microstructure,
)
[14]:
# Use the default material parameters
material_parameters = pulse.HolzapfelOgden.default_parameters()
[15]:
# Select model for active contraction
active_model = "active_strain"
# active_model = "active_stress"
[16]:
# Set the activation
activation = Constant(0.0)
[17]:
# Create material
material = pulse.HolzapfelOgden(
active_model=active_model, parameters=material_parameters, activation=activation
)
[18]:
# Make Dirichlet boundary conditions
def dirichlet_bc(W):
V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
return DirichletBC(V, Constant((0.0, 0.0, 0.0)), fixed)
[19]:
# Make Neumann boundary conditions
neumann_bc = pulse.NeumannBC(traction=Constant(0.0), marker=free_marker)
[20]:
# Collect Boundary Conditions
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))
[21]:
# Create problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
[22]:
# Solve problem
pulse.iterate.iterate(problem, activation, 0.1)
2021-05-12 08:09:22,568 - pulse.iterate - INFO - Iterating....
2021-05-12 08:09:22,569 - pulse.iterate - INFO - Current control: 0.000
2021-05-12 08:09:22,569 - pulse.iterate - INFO - Target: 0.100
[22]:
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 44),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 219),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 255),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 291),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 327),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 359),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 391),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 431),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 463)],
[Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 42),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 217),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 253),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 289),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 325),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 357),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 389),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 429),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 461)])
[23]:
# Get displacement and hydrostatic pressure
u, p = problem.state.split(deepcopy=True)
[24]:
V = dolfin.VectorFunctionSpace(mesh, "CG", 1)
u_int = dolfin.interpolate(u, V)
new_mesh = Mesh(mesh)
dolfin.ALE.move(new_mesh, u_int)
[25]:
from fenics_plotly import plot
fig = plot(mesh, show=False)
fig.add_plot(plot(new_mesh, color="red", show=False))
fig.show()